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Neutrinos in the cosmic ray flux with energies near 1 EeV and above are detectable with the 
Surface Detector array of the Pierre Auger Observatory. We report here on searches through Auger 
data from 1 January 2004 until 20 June 2013. No neutrino candidates were found, yielding a 
limit to the diffuse flux of ultra-high energy neutrinos that challenges the Waxman-Bahcall bound 
predictions. Neutrino identification is attempted using the broad time-structure of the signals 
expected in the SD stations, and is efficiently done for neutrinos of all flavors interacting in the 
atmosphere at large zenith angles, as well as for “Earth-skimming” neutrino interactions in the case 
of tau neutrinos. In this paper the searches for downward-going neutrinos in the zenith angle bins 
60° —75° and 75° —90° as well as for upward-going neutrinos, are combined to give a single limit. The 



4 


90% C.L. single-flavor limit to the diffuse flux of ultra-high energy neutrinos with an i? ^ spectrum 
in the energy range 1.0 x 10^^ eV - 2.5 x 10^® eV is E^dN^/dEi, < 6.4 x 10“® GeV cm“® s“^ sr“^. 

PACS numbers: 95.55.Vj, 95.85.Ry, 98.70.Sa 

Keywords: Ultra-high-energy cosmic rays and neutrinos, high-energy showers, ground detector arrays, Pierre 
Auger Observatory 


I. INTRODUCTION 

The flux of ultra-high energy cosmic rays (UHECRs) 
above ~ 5 x 10^® eV is known to be suppressed with 
respect to that extrapolated from lower energies. This 
feature has been seen in the UHECR spectrum mm, 
with the position of the break being compatible with the 
Greisen-Zatsepin-Kuzmin (GZK) effect [ 3 ], i.e. the in¬ 
teraction of UHECRs with the cosmic microwave back¬ 
ground (CMB) radiation. However, other explanations 
are possible, most prominently a scenario where the lim¬ 
iting energy of the UHECR sources is being observed [ 1 ] . 
Key to distinguishing between these two scenarios is the 
determination of the composition of the UHECRs 
with the second scenario predicting increasing fractions 
of primaries heavier than protons as energy increases @] . 

Above ~ 5 X 10^® eV cosmic-ray protons interact with 
CMB photons and produce ultra-high energy cosmogenic 
neutrinos of energies typically 1/20 of the proton energy 
(3- Their fluxes are uncertain and at EeV energies they 
depend mostly on the evolution with redshift z of the 
unknown sources of UHECRs, and on their spectral fea¬ 
tures at injection. Protons typically produce more neu¬ 
trinos than heavier primaries do mu, so measurement of 
the neutrino flux gives information on the nature of the 
primaries. In this respect the observation of UHE neutri¬ 
nos can provide further hints on the dominant scenario of 
UHECR production [9], as well as on the evolution with 
z of their sources which can help in their identification 

pnni- 

UHE neutrinos are also expected to be produced in the 
decay of charged pions created in the interactions of cos¬ 
mic rays with matter and/or radiation at their potential 
sources, such as Gamma-Ray Bursts or Active Galactic 
Nuclei among others m- In fact, at tens of EeV, neu¬ 
trinos may be the only direct probe of the sources of 
UHECRs at distances farther than ^ 100 Mpc. 

A breakthrough in the field was the recent detection 
with the IceCube experiment of three neutrinos of en¬ 
ergies just above 1 PeV, including a 2 PeV event which 
is the highest-energy neutrino interaction ever observed, 
followed by tens of others above ~ 30 TeV representing 
a ^ 5.7 a excess above atmospheric neutrino background 
|12j . The measured flux is close to the Waxman-Bahcall 


* Now at Fermilab, Batavia, IL, USA 
^ Deceased 

^ Now at CERN, Geneva, Switzerland 
^ Also at Vrije Universiteit Brussels, Belgium 
^ auger_spokespersons@fnal.gov 


upper bound to the UHE neutrino flux [T3j, although 
with a steeper spectrum m- 

In the EeV energy range, i.e. about three orders of 
magnitude above the most energetic neutrinos detected 
in IceCube, neutrinos have so far escaped detection by 
existing experiments. These can be detected with a vari¬ 
ety of techniques m, among them with arrays of particle 
detectors at ground. 

In this work we report on the search for EeV neutrinos 
in data taken with the Surface Detector array (SD) of the 
Pierre Auger Observatory m- A blind scan of data from 
1 January 2004 up to 20 June 2013 has yielded no neu¬ 
trino candidates and an updated and stringent limit to 
the diffuse flux of UHE neutrino flux has been obtained. 


II. SEARCHING FOR UHE NEUTRINOS IN 
AUGER 

The concept for identification of neutrinos is rather 
simple. While protons, heavier nuclei, and even photons 
interact shortly after entering the atmosphere, neutrinos 
can initiate showers quite deep in the atmosphere. At 
large zenith angles the atmosphere is thick enough so 
that the electromagnetic component of nucleonic cosmic 
rays gets absorbed and the shower front at ground level 
is dominated by muons (“old” shower front). On the 
other hand, showers induced by neutrinos deep in the at¬ 
mosphere have a considerable amount of electromagnetic 
component at the ground (“young” shower front). The 
Surface Detector array (SD) of the Pierre Auger Obser¬ 
vatory is not directly sensitive to the muonic and elec¬ 
tromagnetic components of the shower separately, nor to 
the depth at which the shower is initiated. In the ~1600 
water-Cherenkov stations of the SD of the Pierre Auger 
Observatory, spread over an area of ^3000 km^, sepa¬ 
rated by 1.5 km and arranged in a triangular grid, the sig¬ 
nals produced by the passage of shower particles are digi¬ 
tised with Flash Analog to Digital Converters (FADC) 
with 25 ns resolution. This allows us to distinguish nar¬ 
row signals in time induced by inclined showers initiated 
high in the atmosphere, from the broad signals expected 
in inclined showers initiated close to the ground. 

Applying this simple idea, with the SD of the Pierre 
Auger Observatory m we can efficiently detect inclined 
showers and search for two types of neutrino-induced 
showers at energies above about 1 EeV: 

1. Earth-skimming (ES) showers induced by tau neu¬ 
trinos (vr) that travel in a slightly upward direction 
with respect to ground. Vr can skim the Earth’s 
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crust and interact relatively close to the surface in¬ 
ducing a tau lepton which escapes the Earth and 
decays in flight in the atmosphere, close to the SD. 

Typically, only Earth-skimming v-r -induced show¬ 
ers with zenith angles 90° < 9 < 95° may be iden¬ 
tified. 

2. Showers initiated by any neutrino flavor moving 
down at large angles with respect to the vertical 
that interact in the atmosphere close to the sur¬ 
face detector array through charged-current (CC) 
or neutral-current (NC) interactions. We include 
here showers induced by v-r interacting in the moun¬ 
tains surrounding the Pierre Auger Observatory. 
Although this latter process is exactly equivalent 
to the “Earth-skimming” mechanism, it is included 
in this class because such showers are also going 
downwards. In the following we will refer to all 
these types of showers as “downward-going” (DG) 
;/-induced showers. 

With the aid of Monte Carlo simulations we have 
established that this search can be performed effi¬ 
ciently as long as it is restricted to showers with 
zenith angles 6 > 60°. Due to the characteristics 
of these showers depending on the zenith angle, the 
search in this channel was performed in two angu¬ 
lar subranges: (a) “low” zenith angle (DGL) corre¬ 
sponding to 60° < 9 < 75° and (b) “high” zenith 
angle (DGH) with 75° < 6» < 90°. 

A. General procedure 

The identification of potential neutrino-induced show¬ 
ers is based on first selecting those events that arrive in 
rather inclined directions, and then selecting among them 
those with FADC traces that are spread in time, indica¬ 
tive of the early stage of development of the shower and a 
clear signature of a deeply interacting neutrino triggering 
the SD. 

First of all, events occurring during periods of data ac¬ 
quisition instabilities |16j are excluded. For the remain¬ 
ing events the FADC traces of the triggered stations are 
first “cleaned” to remove accidental signals m induced 
mainly by random atmospheric muons arriving closely 
before or after the shower front. These muons are typi¬ 
cally produced in lower energy showers (below the energy 
threshold of the SD of the Auger Observatory) that ar¬ 
rive by chance in coincidence with the triggering shower. 
A procedure to select the stations participating in the 
event described in [13 [18] is then applied, with the event 
accepted if the number of accepted stations Ngt is at least 
three (four) in the Earth-skimming (downward-going) se¬ 
lections. 

From the pattern (footprint) of stations at ground a 
length L along the arrival direction of the event and a 
width W perpendicular to it characterizing the shape of 
the footprint are extracted m- The ratio L/W ~ 1 in 


vertical events, increasing gradually as the zenith angle 
increases. Very inclined events typically have elongated 
patterns on the ground along the direction of arrival and 
hence large values of L(W. A cut in L/W is therefore 
a good discriminator of inclined events. Another indica¬ 
tion of inclined events is given by the apparent speed V of 
the trigger from a station i to a station j, averaged over 
all pairs {i,j) of stations in the event. This observable 
denoted as {V) is obtained from the distance between 
the stations after projection along L and from the differ¬ 
ence in trigger times of the stations. In vertical showers 
{V) exceeds the speed of light since all triggers occur 
at roughly the same time, while in very inclined events 
{V) is concentrated around the speed of light. More¬ 
over its Root-Mean-Square (RMS(V)) value is small. For 
downward-going events only, a cut on the reconstructed 
zenith angle 9rec is applied [l8] . 

Once inclined showers are selected the next step is to 
identify young showers. A Time-over-Threshold (ToT) 
triggeiQ is usually present in SD stations with signals 
extended in time, while narrow signals induce other local 
triggers. Also the Area-over-Peak ratio (AoP), defined 
as the ratio of the integral of the FADC trace to its peak 
value, normalized to the average signal produced by a 
single muon, provides an estimate of the spread-in-time 
of the traces, and serves as an observable to discriminate 
broad from narrow shower fronts. In particular, a cut on 
AoP allows the rejection of background signals induced 
by inclined hadronic showers, in which the muons and 
their electromagnetic products are concentrated within a 
short time interval, exhibiting AoP values close to the one 
measured in signals induced by isolated muons. These 
observables are used by themselves in the search for i/ 
candidates, or combined in a linear Fisher-discriminant 
polynomial depending on the selection as described later 
in this work. 

As a general procedure and to optimize the numeri¬ 
cal values of the cuts and tune the algorithms needed to 
separate neutrino-induced showers from the much larger 
background of hadronic showers, we divided the whole 
data sample (I January 2004 - 20 June 2013) into two 
parts (excluding periods of array instability). A selec¬ 
tion dependent fraction of the data ~ 20%, along with 
Monte Garlo simulations of UHE neutrinos, is dedicated 
to define the selection algorithm, the most efficient ob¬ 
servables and the value of the cuts on them. These data 
are assumed to be overwhelmingly constituted of back¬ 
ground showers. The applied procedure is conservative 
because the presence of neutrinos in the training data 
would result in a more severe definition of the selection 


^ This trigger is intended to select sequences of small signals in 
the FADC traces spread in time. It requires at least 13 bins in 
120 FADC bins of a sliding window of 3 fis above a threshold of 
0.2 (the peak value of the signal expected for a vertical 

muon crossing the station), in coincidence in 2 out of 3 PMTs 

m- 
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criteria. The remaining fraction of data is not used until 
the selection procedure is established, and then it is “un¬ 
blinded” to search for neutrino candidates. We used real 
data to train the selections instead of Monte Carlo sim¬ 
ulations of hadronic showers, the primary reason being 
that the detector simulation may not account for all pos¬ 
sible detector fluctuations that may induce events that 
constitute a background to UHE neutrinos, while they 
are contained in data. It is important to remark that 
this is the same selection procedure and training period 
as in previous publications nmn], which is applied in 
this work to a larger data set. 

Regarding the Monte Carlo simulations, the phase 
space of the neutrino showers reduces to three variables: 
the neutrino energy the incidence zenith angle 9 and 
the interaction depth D in the atmosphere for downward¬ 
going neutrinos, or the altitude he of the r decay above 
ground in the case of Earth-skimming neutrinos. Show¬ 
ers were simulated with energies from log(ET-/eV) = 17 
to 20.5 in steps of 0.5, zenith angles from 90.1° to 95.9° 
in steps of 0.01 rad (ES) and from 60° to 90° in steps of 
0.05 rad (DG). The values of he range from 0 to 2500 m 
(in steps of 100 m) whereas D is uniformly distributed 
along the shower axis in steps of 100 g cm“^. 

We have described the general procedure to search 
for Earth-skimming Vr and downward-going j^-induced 
showers. However the two searches (ES and DG) dif¬ 
fer in several aspects that we describe in the following 
sections. 


B. Earth-skimming (ES) neutrinos 

With Monte Carlo simulations of UHE Vr propagating 
inside the Earth, we have established that r leptons 
above the energy threshold of the SD are efficiently pro¬ 
duced only at zenith angles between 90° and 95°. For 
this reason, in the Earth-skimming analysis we place 
very restrictive cuts to select only quasi-horizontal show¬ 
ers with largely elongated footprints: L/W > 5 and 
{V) S [0.29,0.31] m ns“^ with RMS(U)< 0.08 m ns“^ 
(see Table |l]0 

In the ES selection, the neutrino identification vari¬ 
ables include the fraction of stations with ToT trigger 
and having AoP> 1.4 for data prior to 31 May 2010 [T7]. 
This fraction is required to be above 60% of the trig¬ 
gered stations in the event. The final choice of the val¬ 
ues of these cuts was made by requiring zero background 
events in the training data sample, corresponding to 1% 
of the events recorded up to that date. For data beyond 


^ The axis of Earth-skimming showers travelling in the upward 
direction does not intersect the ground, contrary to the case for 
downward-going showers. For this reason, we exploit the prop¬ 
erties of the footprint generated by the shower particles that 
deviate laterally from the shower axis and trigger the SD water- 
Cherenkov stations. 



Figure 1. Distributions of (AoP) (the variable used to iden- 
tifiy neutrinos in the ES selection for data after 1 June 2010) 
after applying the inclined shower selection in Table 13 Gray- 
filled histogram: the data in the training period. Black his¬ 
togram: data in the search period. These two distributions 
are normalised to the same number of events for compari¬ 
son purposes. Blue histogram: simulated ES u-r events. The 
dashed vertical line represents the cut on (AoP) > 1.83 above 
which a data event is regarded as a neutrino candidate. An 
exponential fit to the tail of the distribution of training data 
is also shown as a red dashed line (see text for explanation). 


1 June 2010 a new methodology and a new set of efficient 
selection criteria was established based on an improved 
and enlarged library of ES simulated Vt events and on 
a larger period of training data. In particular, we used 
the average value of AoP ((AoP)) over all the triggered 
stations in the event as the main observable to discrim¬ 
inate between hadronic showers and ES neutrinos. The 
new methodology allows us to place the value of the cut 
on (AoP) using the tail of its distribution as obtained in 
real data (which was seen to be consistent with an ex¬ 
ponential shape as shown in Fig. [^. This tail was fitted 
and extrapolated to find the value of the cut correspond¬ 
ing to less than 1 expected event per 50 yr on the full 
SD array. As a result, an event is tagged as a neutrino 
candidate if (AoP) > 1.83 (see Tableland Fig. [^. The 
new methodology is not applied to the data prior to 31 
May 2010 since that data period was already unblinded 
to search for UHE neutrinos under the older cuts Hi!. 

Roughly ^ 95% of the simulated inclined v^. events 
producing r leptons above the energy threshold of the SD 
are kept after the cut on (AoP). The search for neutrinos 
is clearly not limited by background in this channel. 


C. Downward-going (DG) neutrinos 

In the high zenith angle range of the downward-going 
analysis (DGH) the values of the cuts to select inclined 
events are obtained in Monte Garlo simulations of events 
with 9 > 75°. Due to the larger angular range compared 
to Earth-skimming less stringent criteria are applied, 
namely L/W > 3, (U) < 0.313 m ns-\ RMS(U)/(U) < 
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Figure 2. Distributions of the Fisher variable T in inclined 
events selected by the “Inclined Showers” DGH criteria in 
Table |T] before applying the “Young Showers” cuts. In par¬ 
ticular the distribution of events with number of triggered 
tanks 7 < N^x, < 11 is shown. Gray-filled histogram: data in 
the training period corresponding to ~ 23% of the whole data 
sample between 1 January 2004 and 20 June 2013. Black line: 
data in the search period. The distributions are normalised 
to the same number of events for comparison purposes. Blue 
line: simulated DGH v events. The dashed vertical line rep¬ 
resents the cut on 7^ > 3.28 above which a data event is 
regarded as a neutrino candidate. The red dashed line repre¬ 
sents an exponential fit to the tail of the training distribution 
(see text for explanation). 


0.08 plus a further requirement that the reconstructed 
zenith angle 0rec > 75° (see [18] and Table |T] for full de¬ 
tails) . 

In the low zenith angle range (DGL) corresponding to 
60° <9 < 75°, LjW, {V) and RMS(y)/(I4) are less effi¬ 
cient in selecting inclined events than the reconstructed 
zenith angle Orec, and for this reason only a cut on ^rec 
is applied, namely 58.5° < 0rec < 76.5°, which includes 
some allowance to account for the resolution in the an¬ 
gular reconstruction of the simulated neutrino events. 

After the inclined shower selection is peformed, the dis¬ 
crimination power is optimized with the aid of the multi¬ 
variate Fisher discriminant method |19] . A linear combi¬ 
nation of observables is constructed which optimizes the 
separation between background hadronic inclined show¬ 
ers occuring during the downward-going training period, 
and Monte Carlo simulated i^-induced showers. The 
method requires as input a set of observables. For that 
purpose we use variables depending on the dimensionless 
Area-over-Peak (AoP) observable - as defined above - of 
the FADC traces. 

In the DGH channel, due to the inclination of the 
shower the electromagnetic component is less attenuated 
at the locations of the stations that are first hit by a deep 
inclined shower [early stations) than in the stations that 
are hit last [late stations). From Monte Carlo simulations 
of induced showers with 9 > 75° we have established 
that in the first few early stations the typical AoP values 
range between 3 and 5, while AoP tends to be closer to 


1 in the late stations. Based on this simple observation 
and as already reported in [18] , we have found a good dis¬ 
crimination when the following ten variables are used to 
construct the linear Fisher discriminant variable IF: the 
AoP and (AoP)^ of the four stations that trigger first in 
each event, the product of the four AoPs, and a global 
parameter that measures the asymmetry between the av¬ 
erage AoP of the early stations and those triggering last 
in the event (see [TS] for further details and Table |l]). 

The selection of neutrino candidates in the zenith an¬ 
gle range 60° < 9 < 75° (DGL) is more challenging since 
the electromagnetic component of background hadronic 
showers at ground increases as the zenith angle decreases 
because the shower crosses less atmosphere before reach¬ 
ing the detector level. Out of all triggered stations of an 
event in this angular range, the ones closest to the shower 
core exhibit the highest discrimination power in terms of 
AoP. In fact it has been observed in Monte Carlo simu¬ 
lations that the first triggered stations can still contain 
some electromagnetic component for background events 
and, for this reason, it is not desirable to use them for dis¬ 
crimination purposes. The last ones, even if they are trig¬ 
gered only by muons from a background hadronic shower, 
can exhibit large values of AoP because they are far from 
the core where muons are known to arrive with a larger 
spread in time. Based on the information from Monte 
Carlo simulations, the variables used in the Fisher dis¬ 
criminant analysis are the individual AoP of the four or 
five stations (depending on the zenith angle) closest to 
the core, and their product [IHl- In the DGL analysis it 
is also required that at least 75% of the triggered stations 
closest to the core have a ToT local trigger [20]. 

Once the Fisher discriminant F is defined, the next 
step is to define a numerical value J^cut that efficiently 
separates neutrino candidates from regular hadronic 
showers. As was done for the variable (AoP) in the 
Earth-skimming analysis, J^cut was fixed using the tail of 
the distribution of F in real data, which is consistent with 
an exponential shape in all cases. An example is shown 
in Fig.[^ The tail was fitted and extrapolated to find the 
value of J^cut corresponding to less than 1 expected event 
per 50 yr on the full SD array [THJ |5D] . Roughly ~ 85% 
(^ 60%) of the simulated inclined v events are kept after 
the cut on the Fisher variable in the DGH (DGL) se¬ 
lections. The smaller efficiencies for the identification of 
neutrinos in the DGL selection are due to the more strin¬ 
gent criteria in the angular bin 9 G (60°, 75°) needed to 
reject the larger contamination from cosmic-ray induced 
showers. 


III. DATA UNBLINDING AND EXPOSURE 
CALCULATION 

A. Data unblinding 

No events survived when the Earth-skimming and 
downward-going selection criteria explained above and 















Selection 

Earth-skimming (ES) 

Downward-going 

Downward-going 



high angle (DGH) 

low angle (DGL) 

Flavours & Interactions 

V-r CC 

CC & NC 

Ve, Vfi , Vt CC & NC 

Angular range 

6 > 90° 

e G (75°,90°) 

9 G (60°, 75°) 

N° of Stations (Mt) 

Mt > 3 

Mt > 4 

Nst > 4 


- 

6>,ec > 75° 

e (58.5°, 76.5°) 

Inclined 

L/W > 5 

L/W > 3 

- 

Showers 

{V) € (0.29, 0.31) m ns’^ 

{V) < 0.313 m ns“^ 

- 


RMS(F) < 0.08 m ns”! 

RMS(F)/(F) < 0.08 

- 


Data: 1 January 2004 - 31 May 2010 


> 75% of stations close to 


> 60% of stations with 


shower core with ToT trigger 

Young 

ToT trigger & AoP >1.4 

Fisher discriminant based 

& 

Showers 

Data: 1 June 2010 - 20 June 2013 

on AoP of early stations 

Fisher discriminant based 


(AoP) > 1.83 


on AoP of early stations 


AoPmin > 1.4 if Aat=3 


close to shower core 


Table I. Observables and numerical values of cuts applied to select inclined and young showers for Earth-skimming and 
downward-going neutrinos. See text for explanation. 


summarized in Table |T] are applied blindly to the data 
collected between 1 January 2004 and 20 June 2013. For 
each selection the corresponding training periods are ex¬ 
cluded from the search. After the unblinding we tested 
the compatibility of the distributions of discriminating 
observables in the search and training samples. Exam¬ 
ples are shown in Fig. for the (AoP) variable in the 
Earth-skimming analysis, and in Fig. for the Fisher 
variable in the DGH analysis. In particular fitting the 
tails of the corresponding distributions to an exponential, 
we obtained compatible parameters within 1 a statistical 
uncertainties. 


B. Exposure calculation 

1. Neutrino identification efficiencies 

The selection criteria in Table |lj were also applied to 
neutrino-induced showers simulated with Monte Carlo, 
and the identification efficiencies cest cdgh, cdgl for 
each channel - defined as the fraction of simulated events 
passing the cuts - were obtained. 

A large set of Monte Carlo simulations of neutrino- 
induced showers was performed for this purpose, cov¬ 
ering the whole parameter space where the efficiency is 
expected to be sizeable. In the case of Earth-skimming 
v-r induced showers, the efficiency depends on the energy 
of the emerging r leptons E^, on the zenith angle 9 and 
on the altitude of the decay point of the r above ground. 
These efficiencies are averaged over azimuthal angle and 
the T decay channels. The maximum efficiency that can 
be reached is 82.6%, the 17.4% remaining corresponds to 
the channel in which the t decays into a /i which is un¬ 


likely to produce a detectable shower close to ground. In 
the case of downward-going neutrinos the identification 
efficiency depends on neutrino flavor, type of interaction 
(CC or NC), neutrino energy zenith angle 9, and dis¬ 
tance D measured from ground along the shower axis at 
which the neutrino is forced to interact in the simula¬ 
tions. 

The identification efficiencies depend also on time, 
through the changing configuration of the SD array that 
was growing steadily since 2004 up to 2008, and be¬ 
cause the fraction of working stations - although typically 
above 95% - is changing continuously with time. Also the 
continuous monitoring of the array reveals a slight evo¬ 
lution with time of the optical properties of the water- 
Cherenkov stations (see below). Although the number 
of working stations and their status are monitored ev¬ 
ery second and as a consequence the SD configuration is 
known with very good accuracy at any instant of time, 
in practice, to avoid having to cope with an impracti- 
cally large number of configurations, different strategies 
were devised to calculate in an accurate and less time- 
consuming manner the actual identification efficiencies 
(as explained in [13 [min])- 

The evolution of the optical properties of the water- 
Cherenkov stations was taken into account in an effective 
way in the calculation of the exposure. The main effect of 
this evolution is a decrease with time of the decay-time 
of the light as obtained from the monitoring data that 
revealed a continuous decrease of ^ 10% from 2004 un¬ 
til the end of the data period used in this work (20 June 
2013). This induces a reduction of the AoP and, as a con¬ 
sequence, the trigger efficiency changes with time. These 
changes were accounted for in the calculation of the ex¬ 
posure by dividing the whole data set into three separate 
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periods and assuming that in each of them the decay-time 
of the light in the tank remained approximately constant 
as seen in data. A conservative approach was adopted 
by choosing constant values of the light decay-time be¬ 
low the actual curve in the three periods. 


8. Combination of selections 

In previous publications HZlIIlllD] the fraction of v- 
induced Monte Carlo events identified as neutrino can¬ 
didates was obtained by applying each particular set of 
selection criteria (ES, DGH, DGL) only to its correspond¬ 
ing set of simulated showers (ES, DGH or DGL). In this 
work the fraction of selected events is further increased 
by applying the three sets of criteria to each sample of 
simulated showers (ES, DGH, DGL) regardless of chan¬ 
nel. With this procedure the fraction of identified Monte 
Garlo events is enhanced as, for instance, an ES simulated 
shower induced by a Vr might not fulfill the requirements 
of the ES selection, but might still pass the DGH or DGL 
criteria, and hence contribute to the fraction of identified 
events. The enhancement in the fraction of events when 
applying this “combined” analysis depends on the par¬ 
ticular set of Monte Carlo simulations. For instance ap¬ 
plying the three criteria to the DGH Monte Carlo sample 
identifies a fraction of neutrino events ~ 1.25 larger than 
when the DGH criteria are applied alone, the enhance¬ 
ment coming mainly from events with 3 stations rejected 
by the DGH criteria but accepted by ES. The applica¬ 
tion of the three criteria to the ES Monte Carlo sample 
however results in a smaller enhancement ~ 1.04. 


3. Exposure calculation 

For downward-going neutrinos, once the efficiencies 
euG{Ev,S, D,t) are obtained, the calculation of the ex¬ 
posure involves folding them with the SD array aperture 
and the v interaction probability at a depth D for a neu¬ 
trino energy This calculation also includes the pos¬ 
sibility that downward-going Vr interact with the moun¬ 
tains surrounding the Observatory. Integrating over the 
parameter space except for Ei, and in time over the search 
periods and summing over all the interaction channels 
yields the exposure [iHiiin]. 

In the Earth-skimming channel, eEs(£'r, 0, Xu) are also 
folded with the aperture, with the probability density 
function of a tau emerging from the Earth with energy E^. 
(given a neutrino with energy E^, crossing an amount of 
Earth determined by the zenith angle 9), as well as with 
the probability that the r decays at an altitude h, m- 
An integration over the whole parameter space except for 
El, and time gives the exposure [17] . 

The exposures £eS) ^dgh and Edgl obtained for the 
search periods of each selection are plotted in Fig.j^along 
with their sum £tot- The exposure to Earth-skimming 
neutrinos is higher than that to downward-going neu¬ 


trinos, partially due to the longer search period in the 
Earth-skimming analysis, and partially due to the much 
larger neutrino conversion probability in the denser tar¬ 
get of the Earth’s crust compared to the atmosphere. The 
larger number of neutrino flavors and interaction chan¬ 
nels that can be identified in the DGH and DGL analysis, 
as well as the broader angular range 60° < 9 < 90° partly 
compensates the dominance of the ES channel. The ES 
exposure flattens and then falls above ^ 10^® eV as there 
is an increasing probability that the r decays high in the 
atmosphere producing a shower not triggering the array, 
or even that the r escapes the atmosphere before decay¬ 
ing. At the highest energies the DGH exposure domi¬ 
nates. The DGL exposure is the smallest of the three, 
mainly due to the more stringent criteria needed to ap¬ 
ply to get rid of the larger background nucleonic showers 
in the zenith angle bin 60° < 9 < 75°. 

The relative contributions of the three channels to 
the total expected event rate for a differential flux 
behaving with energy as dNy{Ey)/dE,, oc are 

ES:DGH:DGL~0.84:0.14:0.02 respectively, where the 
event rate is obtained as: 

iVevt= / ^{Ei,) dE„ (1) 

Je„ dEu 

C. Systematic uncertainties 

Several sources of systematic uncertainty have been 
considered. Some of them are directly related to the 
Monte Garlo simulation of the showers, i.e., generator 
of the neutrino interaction either in the Earth or in the 
atmosphere, parton distribution function, air shower de¬ 
velopment, and hadronic model. 

Other uncertainties have to do with the limitations on 
the theoretical models needed to obtain the interaction 
cross-section or the r energy loss at high energies. In the 
Earth-skimming analysis the model of energy loss for the 
T is the dominant source of uncertainty, since it deter¬ 
mines the energy of the emerging rs after propagation 
in the Earth; the impact of this on the downward-going 
analysis is much smaller since r energy losses are only 
relevant for Vr interacting in the mountains, a channel 
that is estimated to contribute only ^ 15% to the DGH 
exposure HHj. 

The uncertainty on the shower simulation, that stems 
mainly from the different shower propagation codes and 
hadronic interaction models that can be used to model 
the high energy collisions in the shower, contributes sig¬ 
nificantly in the ES and DG channels. 

The presence of mountains around the Observatory - 
which would increase the target for neutrino interactions 
in both cases - is explicitly simulated and accounted for 
when obtaining the exposure of the SD to downward¬ 
going neutrino-induced showers, and as a consequence 
does not contribute directly to the systematic uncertain¬ 
ties. However, it is not accounted for in the Earth- 
skimming channel and instead we take the topography 
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around the Observatory as a source of systematic uncer¬ 
tainty. 

In the three channels the procedure to incorporate the 
systematic uncertainties is the same. Different combi¬ 
nations of the various sources of systematic uncertainty 
render different values of the exposure and a systematic 
uncertainty band of relative deviation from a reference 
exposure (see below) can be constructed for each chan¬ 
nel and for each source of systematic uncertainty. For 
a given source of uncertainty the edges of the ES, DGH 
and DGL bands are weighted by the relative importance 
of each channel as given before and added linearly or 
quadratically depending on the source of uncertainty. In 
Table [T^ we give the dominant sources of systematic un¬ 
certainty and their corresponding combined uncertainty 
bands obtained in this way. The combined uncertainty 
band is then incorporated in the value of the limit itself 
through a semi-Bayesian extension of the Feldman- 
Cousins approach [22] . 

In the calculation of the reference exposure the v- 
nucleon interaction in the atmosphere for DG neutrinos 
(including GG and NG channels) is simulated with HER- 
WIG [3S]. In the case of i^r CC interactions, a dedicated, 
fast and flexible code is used to simulate the r lepton 
propagation in the Earth and/or in the atmosphere. The 
T decay is performed with the TAUOLA package |M|. In 
all cases we adopted the nucleon cross-section in [25] . 
In a second step, the AIRES code [26] is used to simulate 
the propagation in the atmosphere of the particles pro¬ 
duced in the high energy v interaction or in the r lepton 
decay. The types, energies, momenta and times of the 
particles reaching the SD level are obtained. The last 
stage is the simulation of the SD response (PMT signals 
and FADG traces). This involves a modification of the 
“standard” sampling procedure in m to regenerate par¬ 
ticles in the SD stations from the “thinned” air shower 
simulation output, that was tailored to the highly in¬ 
clined showers involved in the search for neutrinos. Light 
production and propagation inside the station is based on 
GEANT4 [3H] with the modifications to account for the 
evolution of the light decay-time explained above. These 
two latter changes roughly compensate each other, with 
the net result being a few percent decrease of the ex¬ 
posure with respect to that obtained with the standard 
thinning procedure and a constant average value of the 
light decay-time. 


IV. RESULTS 

Using the combined exposure in Fig. [^ and assuming a 
differential neutrino flux dN{E^)/dE^ = k ■ as well 
as a Vg : : Vt = 1 '■ I '■ 1 flavor ratio, an upper limit on 

the value of k can be obtained as: 


/e„ ^tot(E'y) dE^, 


Source of systematic 

Combined uncertainty band 

Simulations 

~ +4%, -3% 

n cross section & r E-loss 

~ -b34%, -28% 

Topography 

~ -bl5%, 0% 

Total 

~ -b37%, -28% 


Table II. Main sources of systematic uncertainties and their 
corresponding combined uncertainty bands (see text for de¬ 
tails) representing the effect on the event rate defined in 
Eq. Q. The uncertainty due to “Simulations” includes: inter¬ 
action generator, shower simulation, hadronic model, thinning 
and detector simulator. The uncertainty due to “r energy- 
loss” affects the ES channel and also the DGH but only to 
Vt with 0 > 88° going through the mountains surrounding 
the Pierre Auger Observatory. However it does not affect the 
DGL channel. The topography around the Observatory is not 
accounted for in the ES channel and is taken as a systematic 
uncertainty that would increase the event rate. 



Figure 3. Combined exposure of the SD of the Pierre Auger 
Observatory (1 January 2004 - 20 June 2013) as a function 
of neutrino energy after applying the three sets of selection 
criteria in Table [I] to Monte Carlo simulations of UHE neu¬ 
trinos (see text for explanation). Also shown are the individ¬ 
ual exposures corresponding to each of the three selections. 
For the downward-going channels the exposure represents the 
sum over the three neutrino flavors as well as CC and NC 
interactions. For the Earth-Skimming channel, only Vt CC 
interactions are relevant. 


The actual value of the upper limit on the signal events 
(A^up) depends on the number of observed events (0 in 
our case) and expected background events (conserva¬ 
tively assumed to be 0), as well as on the confidence 
level required (90% G.L. in the following). Using a semi- 
Bayesian extension m of the Feldman-Cousins approach 
[22] to include the uncertainties in the exposure we ob- 
tair0 fVup = 2.39. The single-flavor 90% G.L. limit is: 


^ To calculate A^up we use POLEH—h m- The signal efficiency un¬ 
certainty is 0.19 with an asymmetric band (see Table [Till. This 
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kgo < 6.4 X 10 ® GeV cm ^ s ^ sr (3) 

The limit applies in the energy interval ^ 1.0 x 10^^ eV — 
2.5 X 10^® eV where the cumulative number of events as 
a function of neutrino energy increases from 5% to 95% 
of the total number, i.e. where ~ 90% of the total event 
rate is expected. It is important to remark that this 
is the most stringent limit obtained so far with Auger 
data, and it represents a single limit combining the three 
channels where we have searched for UHE neutrinos. The 
limit to the flux normalization in Eq. (§ is obtained in¬ 
tegrating the denominator of Eq. (§ in the whole energy 
range where Auger is sensitive to UHE neutrinos. This 
is shown in Fig. , along with the 90% C.L. limits from 
other experiments as well as several models of neutrino 
flux production (see caption for references). The denom¬ 
inator of Eq. ([^ can also be integrated in bins of energy, 
and a limit on k can also be obtained in each energy bin 
|35| . This is displayed in Fig. where the energy bins 
have a width of 0.5 in log^o and where we also show 
the whole energy range where there is sensitivity to neu¬ 
trinos. The limit as displayed in Fig. [fallows us to show 
at which energies the sensitivity of the SD of the Pierre 
Auger Observatory peaks. 

The search period corresponds to an equivalent of 6.4 
years of a complete Auger SD array working continuously. 
The inclusion of the data from 1 June 2010 until 20 June 
2013 in the search represents an increase of a factor ~ 1.8 
in total time quantified in terms of equivalent full Auger 
years with respect to previous searches [IIlIIH]. Further 
improvements in the limit come from the combination of 
the three analysis into a single one, using the procedure 
explained before that enhances the fraction of identified 
neutrinos especially in the DGH channel. 

In Table [m] we give the expected total event rates for 
several models of neutrino flux production. 

Several important conclusions and remarks can be 
stated after inspecting Figs. and [^and Table [In| 

1. The maximum sensitivity of the SD of the 
Auger Observatory is achieved at neutrino energies 
around EeV, where most cosmogenic models of v 
production also peak (in a x dN/dE^ plot). 

2. The current Auger limit is a factor ~ 4 below the 
Waxman-Bahcall bound on neutrino production in 
optically thin sources [13]. The SD of the Auger 
Observatory is the first air shower array to reach 
that level of sensitivity. 

3. Some models of neutrino production in astrophys- 
ical sources such as Active Galactic Nuclei (AGN) 
are excluded at more than 90% G.L. For the model 


yields a value of A^up = 2.39 slightly smaller than the nominal 
2.44 of the Feldman-Cousins approach. 
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Figure 4. Top panel: Upper limit (at 90% C.L.) to the 
normalization of the diffuse flux of UHE neutrinos as given 
in Eqs. ([^ and (§, from the Pierre Auger Observatory. We 
also show the corresponding limits from ANITAII [5^ and 
IceCube |.30| experiments, along with expected fluxes for sev¬ 
eral cosmogenic neutrino models that assume pure protons 
as primaries |311133| as well as the Waxman-Bahcall bound 
m- All limits and fluxes converted to single flavor. We used 
Nup = 2.39 in Eq. § to obtain the limit (see text for de¬ 
tails). Bottom panel: Same as top panel, but showing several 
cosmogenic neutrino models that assume heavier nuclei as pri¬ 
maries, either pure iron m or mixed primary compositions 
0 . 


^2 shown in Fig. 14 of [32] we expect ^7 neutrino 
events while none was observed. 

4. Cosmogenic v models that assume a pure primary 
proton composition injected at the sources and 
strong (FRII-type) evolution of the sources are 
strongly disfavored by Auger data. An example 
is the upper line of the shaded band in Fig. 17 in 
[31] (also depicted in Figs. and [^, for which ^4 
events are expected and as consequence that flux 
is excluded at ~98% C.L. Models that assume a 
pure primary proton composition and use the GeV 
7 -ray flux observations by the Fermi-LAT satellite 
detector as an additional constraint, are also dis¬ 
favored. For instance for the model shown as a 
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Diffuse flux 

Expected number of events 

Probability of 

Neutrino Model 

(1 January 2004 - 20 June 2013) 

observing 0 

Cosmogenic - proton, FRII |31| 

~ 4.0 

~ 1.8 X 10"^ 

Cosmogenic - proton, SFR [ST] 

~ 0.9 

~ 0.4 

Cosmogenic - proton, Fermi-LAT, Emin = IQi® eV [33] 

~ 3.2 

~ 4 X 10"^ 

Cosmogenic - proton, Fermi-LAT, Emin = eV [33] 

~ 1.6 

~ 0.2 

Cosmogenic - proton or mixed, SFR & CRB [3] 

~ 0.5 - 1.4 

~0.6 - 0.2 

Cosmogenic - iron, FRII [31] 

~ 0.3 

~ 0.7 

Astrophysical u (AGN) [33] 

~ 7.2 

~ 7 X 10-1 

Exotic [33] 

~ 31.5 

~ 2 X 10-11 


Table III. Number of expected events Nevt in Eq. Q for several theoretical models of UHE neutrino production (see Figs. 
and|^, given the combined exposure of the surface detector array of the Pierre Auger Observatory plotted in Fig. The last 
column gives the Poisson probability exp(—Aevt) of observing 0 events when the number of expected events is Nevt given in 
the second column. 



Figure 5. Upper limit to the normalization of the diffuse 
flux of UHE neutrinos (at 90% C.L. and in bins of width 0.5 
in logj^Q - see text for details) from the Pierre Auger Ob¬ 
servatory (straight steps). We also show the corresponding 
limits from ANITAII (dot-dashed line) and IceCube |30| 
(dashed line) experiments (with appropriate normalizations 
to take into account the energy bin width, and to convert 
to single flavor), along with expected fluxes for several cos- 
mogenic neutrino models [9l 1311133| as well as the Waxman- 
Bahcall bound [13] (all converted to single flavor). 


solid line in the bottom right panel of Fig. 5 in [ 33 ] 
(also depicted in Figs. |^and[^in this work), cor¬ 
responding to the best-fit to the cosmic-ray spec¬ 
trum as measured by HiRes, we expect ^3.2 events. 
As a consequence that model is excluded at more 
than 90% C.L. For this particular model we also 
show in Figs, [d] and the 99% C.L. band result¬ 
ing from the fitting to the HiRes spectrum down to 
L^min = 10^® eV. The Auger limit is also approach¬ 
ing the solid line in the upper left panel of Fig. 5 
in [33j , a model that assumes extragalactic protons 
above i?min = eV [ 53 ], for which ~ 1.6 events 


are expected (see Table III). The Auger direct lim¬ 
its on cosmogenic neutrinos are also constraining 
part of the region indirectly bounded by Fermi-LAT 
observations. 


5. The current Auger limit is less restrictive with 
the cosmogenic neutrino models represented by the 
gray shaded area in the bottom panel of Fig. 
(~ 0.5 t o ~1.4 events are expected as shown in Ta¬ 
ble III) which brackets the lower fluxes predicted 
under a range of assumptions for the composition 
of the primary flux (protons or mixed), source evo¬ 
lution and model for the transition from Galactic to 
extragalactic cosmic-rays The same remark ap¬ 
plies to models that assume pure-iron composition 
at the sources. A 10-fold increase in the current ex¬ 
posure will be needed to reach the most optimistic 
predictions of cosmogenic neutrino fluxes if the pri¬ 
maries are pure iron, clearly out of the range of the 
current configuration of the Auger Observatory. 


6 . A large range of exotic models of neutrino produc¬ 
tion [34] are excluded with C.L. larger than 99%. 

7. In IceCube, neutrino fluxes in the 30 TeV to 2 PeV 
energy range have shown a ~5.7(T excess compared 
to predicted atmospheric neutrino fluxes | 12 j . A 
refinement of the IceCube search technique to ex¬ 
tend the neutrino sensitivity down to 10 TeV [37] . 
yielded a power-law fit to the measured flux with¬ 
out cut-off given by dN/dE = ^o{E,j/Eo)~'^ with 
$0 = 2.06x10-1® GeV-i cm-^ s"! sr"!, Eq = 10® 
GeV, and 7 = 2.46. If this flux is extrapolated to 
10^° eV it would produce ^0.1 events in Auger. 


V. ACKNOWLEDGMENTS 

The successful installation, commissioning, and oper¬ 
ation of the Pierre Auger Observatory would not have 
been possible without the strong commitment and effort 























13 


from the technical and administrative staff in Malargiie. 
We are very grateful to the following agencies and orga¬ 
nizations for financial support: 

Comision Nacional de Energia Atomica, Fundacion 
Antorchas, Gobierno de la Provincia de Mendoza, Mu- 
nicipalidad de Malargiie, NDM Holdings and Valle Las 
Lehas, in gratitude for their continuing cooperation over 
land access, Argentina; the Australian Research Coun¬ 
cil; Conselho Nacional de Desenvolvimento Cientihco e 
Tecnologico (CNPq), Financiadora de Estudos e Pro- 
jetos (FINEP), Fundagao de Amparo a Pesquisa do 
Estado de Rio de Janeiro (FAPERJ), Sao Paulo Re¬ 
search Foundation (FAPESP) Grants No. 2010/07359- 
6 and No. 1999/05404-3, Ministerio de Ciencia e 
Tecnologia (MCT), Brazil; Grant No. MSMT-CR 
LG13007, No. 7AMB14AR005, and the Czech Sci¬ 
ence Foundation Grant No. 14-17501S, Czech Re¬ 
public; Centre de Calcul IN2P3/CNRS, Centre Na¬ 
tional de la Recherche Scientifique (CNRS), Conseil 
Regional Ile-de-France, Departement Physique Nucleaire 
et Corpusculaire (PNC-IN2P3/CNRS), Departement 
Sciences de I’Univers (SDU-INSU/CNRS), Institut La¬ 
grange de Paris (ILP) Grant No. LABEX ANR-10- 
LABX-63, within the Investissements d’Avenir Pro¬ 
gramme Grant No. ANR-ll-IDEX-0004-02, France; 
Bundesministerium fiir Bildung und Forschung (BMBF), 
Deutsche Forschungsgemeinschaft (DFG), Finanzmin- 
isterium Baden-Wiirttemberg, Helmholtz Alliance for 
Astroparticle Physics (HAP), Helmholtz-Gemeinschaft 
Deutscher Forschungszentren (HGF), Ministerium fiir 
Wissenschaft und Forschung, Nordrhein Westfalen, Min¬ 
isterium fiir Wissenschaft, Forschung und Kunst, Baden- 
Wiirttemberg, Germany; Istituto Nazionale di Fisica Nu- 
cleare (INFN), Ministero dell’Istruzione, dell’Universita 
e della Ricerca (MIUR), Gran Sasso Center for As¬ 
troparticle Physics (CFA), CETEMPS Center of Ex¬ 


cellence, Ministero degli Affari Esteri (MAE), Italy; 
Consejo Nacional de Ciencia y Tecnologia (CONA- 
CYT), Mexico; Ministerie van Onderwijs, Cultuur 
en Wetenschap, Nederlandse Organisatie voor Weten- 
schappelijk Onderzoek (NWO), Stichting voor Funda- 
menteel Onderzoek der Materie (FOM), Netherlands; 
National Centre for Research and Development, Grants 
No. ERA-NET-ASPERA/Ol/lI and No. ERA-NET- 
ASPERA/02/I1, National Science Centre, Grants No. 
20I3/08/M/ST9/00322, No. 20I3/08/M/ST9/00728 
and No. HARMONIA 5 - 2013/10/M/ST9/00062, 
Poland; Portuguese national funds and FEDER funds 
within Programa Operacional Factores de Competitivi- 
dade through Fundagao para a Ciencia e a Tecnologia 
(COMPETE), Portugal; Romanian Authority for Sci¬ 
entific Research ANCS, CNDI-UEFISCDI partnership 
projects Grants No. 20/2012 and No. 194/2012, Grants 
No. 1/ASPERA2/2012 ERA-NET, No. PN-H-RU- 
PD-2011-3-0145-17 and No. PN-H-RU-PD-2011-3-0062, 
the Minister of National Education, Programme Space 
Technology and Advanced Research (STAR), Grant No. 
83/2013, Romania; Slovenian Research Agency, Slove¬ 
nia; Comunidad de Madrid, FEDER funds, Ministerio 
de Educacion y Ciencia, Xunta de Galicia, European 
Community 7th Framework Program, Grant No. FP7- 
PEOPLE-2012-IEF-328826, Spain; Science and Tech¬ 
nology Facilities Council, United Kingdom; Depart¬ 
ment of Energy, Contracts No. DE-AC02-07CH11359, 
No. DE-FR02-04ER41300, No. DE-FG02-99ER41107 
and No. de-sc0011689, National Science Foundation, 
Grant No. 0450696, The Grainger Foundation, USA; 
NAFOSTED, Vietnam; Marie Curie-IRSES/EPLANET, 
European Particle Physics Latin American Network, 
European Union 7th Framework Program, Grant No. 
PIRSES-2009-GA-246806; and UNESCO. 


[1] R. U. Abbasi et al. [HiRes Collaboration], Phys. Rev. 
Lett. 100, 101101 (2008). 

[2] J. Abraham et al. [Pierre Auger Collaboration], Phys. 
Rev. Lett. 101, 061101 (2008); J. Abraham et al. [Pierre 
Auger Collaboration], Phys. Lett. B 685, 239 (2010). 

[3] K. Greisen, Phys. Rev. Lett. 16, 748 (1966); 

G. T. Zatsepin and V. A. Kuzmin, Pis’ma ZhZhurnal 
Eksperimental noi i Teoreticheskoi Fiziki 4, 114 (1966), 
JETP Lett. 4, 78 (1966) (Engl. TransL). 

[4] D. Allard, Astropart. Phys. 33, 39 (2012). 

[5] P. Abreu et al. [Pierre Auger Collaboration], Phys. Rev. 
Lett. 104, 091101 (2010); A. Aab et al. [Pierre Auger 
Collaboration], arXiv: 1409.4809 [astro-ph]. 

[6] R. U. Abbassi et al. [Telescope Array Collaboration] As¬ 
tropart. Phys. 64, 49 (2015). 

[7] V. S. Berezinsky and G. T. Zatsepin, Phys. Lett. B 28, 
423 (1969). 

[8] D. Hooper et al, Astropart. Phys. 23, 11 (2005) 11; M. 
Ave et al., Astropart. Phys. 23, 19 (2005) 19. 

[9] D. Allard et al, JCAP 10, 013 (2010). 


[10] D. Seckel and T. Stanev, Phys. Rev. Lett. 95, 141101 
(2005). 

[11] J. K. Becker, Phys. Rep. 458, 173 (2008). 

[12] M. G. Aartsen et al. [IceGube Collab.] Phys. Rev. Lett. 
113, 101101 (2014). 

[13] E. Waxman and J. N. Bahcall, Phys. Rev. D 64, 023002 

( 2001 ). 

[14] B. Baret and V. Van Elewyck, Rep. Prog. Phys. 74, 
046902 (2011). 

[15] J. Abraham et al. The Pierre Auger Collab., Nucl. lu¬ 
strum. Meth. A 523, 50 (2004). 

[16] J. Abraham et al. [Pierre Auger Collaboration], Nucl. 
Instrum. Meth. A 613, 29 (2010). 

[17] J. Abraham et al. [Pierre Auger Collab.], Phys. Rev. Lett. 
100, 211101 (2008); J. Abraham et al. [Pierre Auger 
Collab.], Phys. Rev. D 79, 102001 (2009); P. Abreu et 
al. [Pierre Auger Collab.], Astrophys. J. Lett. 755, L4 
( 2012 ). 

[18] P. Abreu et al. [Pierre Auger Collab.], Phys. Rev. D 84, 
122005 (2011). 



14 


[19] R. Fisher, Ann. Eugenics 7, 179 (1936). 

[20] J. L. Navarro, PhD Thesis, Univ. Granada, Spain (2012). 

[21] J. Conrad et al, Phys. Rev. D 67, 012002 (2003). 

[22] G. J. Feldman, R. D. Cousins, Phys. Rev. D 57, 3873 
(1998). 

[23] G. Corcella et al, HERWIG 6.5, JHEP 01, 010 (2001). 

[24] S. Jadach et al., Comput. Phys. Common. 76, 361 (1993). 

[25] A. Cooper-Sarkar, S. Sarkar, JHEP 0801, 075 (2008). 

[26] S. J. Sciutto, arXiv:astro-ph/9911331 

http://www2.fisica.unlp.edu.ar/auger/aires/ppal.html 

[27] P. Billoir, Astropart. Phys. 30, 270 (2008). 

[28] http://geant4.cern.ch/ 

[29] P. W. Gorham et al. [ANITA Collab.], Phys. Rev. D 85, 
049901(E) (2012). 


[30] M. G. Aartsen et al. [IceCube Collab.] Phys. Rev. D 88, 
112008 (2013). 

[31] K. -H. Kampert, M. Unger, Astropart. Phys. 35, 660 

( 2012 ). 

[32] J. K. Becker, P.L. Biermann, W. Rhode, Astropart. Phys. 
23, 355 (2005). 

[33] M. Ahlers et al., Astropart. Phys. 34, 106 (2010). 

[34] O. Kalashev et al., Phys. Rev. D 66, 063004 (2002). 

[35] L. A. Anchordoqui et al., Phys. Rev. D 66, 103002 (2002) 

[36] V. Berezinsky, A. Z. Gazizov, S. I. Grigorieva, Phys. Rev. 
D 74, 043005 (2006). 

[37] M. G. Aartsen et al. [IceCube Collab.] Phys. Rev. D 91, 
022001 (2015). 



